##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(Seurat)

##########################################################################################
option_list <- list(
    make_option(c("--geneset_file"), type = "character"),
    make_option(c("--magic_exp_file"), type = "character"),
    make_option(c("--rna_data_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){

    ## 表达文件
    rna_data_file <- "~/20231121_singleMuti/results/qc_atac_v3/germ/testis_combined.annotationCellType.qc.Rdata"

    ## 
    geneset_file <- "~/20231121_singleMuti/config/Human_reported_TF2.addInfo.csv"

    ## 填补以后的中位表达文件
    magic_exp_file <- "~/20231121_singleMuti/results/qc_atac_v3/germ/GeneExpression.MeanByCellType.magic.tsv"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/qc_atac_v3/germ/"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

geneset_file <- opt$geneset_file
magic_exp_file <- opt$magic_exp_file
rna_data_file <- opt$rna_data_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################
## 导入数据
a <- load(rna_data_file)
dat_geneset <- data.frame(fread(geneset_file , header = T))
dat_magic <- data.frame(fread(magic_exp_file , header = T))

###########################################################################################
## 计算表达大于0的比例

exp_numeric <- apply(scrnat@assays$RNA@data , 1 , as.numeric)
rownames(exp_numeric) <- colnames(scrnat@assays$RNA@data)
exp_numeric <- t(exp_numeric)

## 所有细胞
all_ratio <- apply(exp_numeric , 1 , function(x){length(which(x>0))})
all_ratio <- data.frame( all = all_ratio/ncol(exp_numeric) )

## 不同类型的细胞
sco_ratio <- sapply(unique(scrnat$cell_type),function(x){
    print(x)
    sapply(unique(rownames(scrnat)),function(y){
        length(which(as.numeric(as.vector(exp_numeric[y,which(scrnat$cell_type==x)])) > 0))/length(which(scrnat$cell_type==x))
    })
})

## 合并
result <- data.frame(cbind( all_ratio , sco_ratio ))

###########################################################################################
## 标记该基因是否在任一细胞满足阈值标注
for( ratio_t in seq(0.01 , 0.25 , 0.01) ){
  print(ratio_t)
  ratio_use <- which(apply( result[,1:10] , 1 , function(x){max(x) >= ratio_t} ))

  result$tmp <- ifelse( rownames(result) %in% names(ratio_use) , "TRUE" , "FALSE" )
  colnames(result)[ncol(result)] <- ratio_t
}

result$gene <- rownames(result)
result <- merge(result , dat_geneset , by.x = "gene" , by.y = "Human_reported_TF" , all.x = T )

out_file <- paste0(out_path , "/GeneExpression.ratio.tsv")
write.table( result , out_file , row.names = F , sep = "\t" , quote = F )

###########################################################################################
## 标记填补以后所在细胞的最大表达值
dat_magic$max_exp <- apply( dat_magic[,c(2:10)] , 1 , max )
for( ratio_t in seq(0.05 , 0.5 , 0.05) ){
  print(ratio_t)
  dat_magic$tmp <- ifelse( dat_magic$max_exp >= ratio_t , "TRUE" , "FALSE" )
  colnames(dat_magic)[ncol(dat_magic)] <- ratio_t
}

dat_magic <- merge(dat_magic , dat_geneset , by.x = "gene" , by.y = "Human_reported_TF" , all.x = T )

out_file <- paste0(out_path , "/GeneExpression.magic.max_mean.tsv")
write.table( dat_magic , out_file , row.names = F , sep = "\t" , quote = F )
